This notebook outlines methods to calculate barometric efficiency from groundwater well data using Python. First we import some sample data, then we analyze the data with a few different techniques.
In [2]:
%matplotlib inline
import pandas as pd
import numpy as np
import os
import glob
import urllib2
import xmltodict
import matplotlib.pyplot as plt
import statsmodels.tsa.tsatools as tools
from pandas.stats.api import ols
from pylab import rcParams
rcParams['figure.figsize'] = 15, 10
In [3]:
def readxle(infile):
def getfilename(path):
# this function extracts the file name without file path or extension
return path.split('\\').pop().split('/').pop().rsplit('.', 1)[0]
# accomodate for files online
if infile[0:3]=='htt' or infile[0:3]=='ftp':
response = urllib2.urlopen(infile)
data = response.read()
obj = xmltodict.parse(data,encoding="ISO-8859-1")
else:
# open text file
with open(infile) as fd:
# parse xml
obj = xmltodict.parse(fd.read(),encoding="ISO-8859-1")
# navigate through xml to the data
wellrawdata = obj['Body_xle']['Data']['Log']
# convert xml data to pandas dataframe
f = pd.DataFrame(wellrawdata)
# get header names and apply to the pandas dataframe
f[obj['Body_xle']['Ch1_data_header']['Identification']] = f['ch1']
f[obj['Body_xle']['Ch2_data_header']['Identification']] = f['ch2']
# add extension-free file name to dataframe
f['name'] = getfilename(infile)
# combine Date and Time fields into one field
f['DateTime'] = pd.to_datetime(f.apply(lambda x: x['Date'] + ' ' + x['Time'], 1))
f[obj['Body_xle']['Ch1_data_header']['Identification']] = f[obj['Body_xle']['Ch1_data_header']['Identification']].convert_objects(convert_numeric=True)
f = f.reset_index()
f = f.set_index('DateTime')
f = f.drop(['Date','Time','@id','ch1','ch2','index','ms'],axis=1)
return f
In [4]:
# import relevant data from github
bpfileurl = "https://raw.githubusercontent.com/inkenbrandt/Earth_Tides/master/pw03%20baro%202015-03-04.xle"
wlfileurl = "https://raw.githubusercontent.com/inkenbrandt/Earth_Tides/master/ag13a%202015-03-03.xle"
bp = readxle(bpfileurl)
bp['bp'] = bp['Level']
wl = readxle(wlfileurl)
# unvented transducer provides absolute pressure in feet of water above transducer
wl['wl_and_bp'] = wl['LEVEL']
# merge bp and wl data
data = pd.merge(bp, wl, left_index=True, right_index=True, how='inner')
# calculate vented water level in feet of water above transducer
data['wl'] = data['wl_and_bp']-data['bp']
# drop old fields
data.drop(['Level','LEVEL','Temperature','TEMPERATURE','name_x','name_y','wl_and_bp'], axis=1, inplace=True)
data.dropna(inplace=True)
In [5]:
# plot the data
x = data.index.to_datetime()
y1 = data['wl']
y2 = data['bp']
fig, ax1 = plt.subplots()
plt.title('Barometric Pressure and Water Level over Time')
ax2 = ax1.twinx()
ax1.plot(x,y1,c='r',label='Water Level')
ax1.set_ylabel('Well Transducer Pressure (ft)', color='r')
ax2.set_ylabel('Barometric Pressure (ft)', color='b')
ax2.plot(x,y2,c='b',label='Barometric Pressure')
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, loc=3)
Out[5]:
This notebook outlines methods to calculate barometric efficiency from groundwater well data. Barometric efficiency is the ratio of changes in well water levels to changes in barometric pressure (Rasmussen and Crawford, 1997). $$ be = \frac{\Delta wl}{\Delta bp} $$
There are several different methods that one can apply to determine barometric efficiency. The most simple is to take first differences of barometric pressure and water level and then compute a linear regression of the two.
In [6]:
data['dwl'] = data['wl'].diff()
data['dbp'] = data['bp'].diff()
regression = ols(y=data['dwl'], x=data['dbp'])
m = regression.beta.x
b = regression.beta.intercept
r = regression.r2
#r = (regression.beta.r)**2
plt.scatter(y=data['dwl'], x=data['dbp'])
y_reg = [data['dbp'][i]*m+b for i in range(len(data['dbp']))]
plt.plot(data['dbp'],y_reg,
label='Regression: Y = {m:.4f}X + {b:.5}\nr^2 = {r:.4f}\n BE = {be:.2f} '.format(m=m,b=b,r=r,be=m))
plt.legend()
plt.xlabel('Sum of Barometric Pressure Changes (ft)')
plt.ylabel('Sum of Water-Level Changes (ft)')
Out[6]:
Clark's (1967) method is one way to calculate barometric efficiency. Merritt (2004) outlines how to perform this method very well.
In [7]:
# clark's method
def clarks(data,bp,wl):
data['dwl'] = data[wl].diff()
data['dbp'] = data[bp].diff()
data['beta'] = data['dbp']*data['dwl']
data['Sbp'] = np.abs(data['dbp']).cumsum()
data['Swl'] = data[['dwl','beta']].apply(lambda x: -1*np.abs(x[0]) if x[1]>0 else np.abs(x[0]), axis=1).cumsum()
plt.figure()
plt.plot(data['Sbp'],data['Swl'])
regression = ols(y=data['Swl'], x=data['Sbp'])
m = regression.beta.x
b = regression.beta.intercept
r = regression.r2
y_reg = [data.ix[i,'Sbp']*m+b for i in range(len(data['Sbp']))]
plt.plot(data['Sbp'],y_reg,
label='Regression: Y = {m:.4f}X + {b:.5}\nr^2 = {r:.4f}\n BE = {be:.2f} '.format(m=m,b=b,r=r,be=m))
plt.legend()
plt.xlabel('Sum of Barometric Pressure Changes (ft)')
plt.ylabel('Sum of Water-Level Changes (ft)')
data.drop(['dwl','dbp','Sbp','Swl'],axis=1,inplace=True)
return m,b,r
m,b,r = clarks(data, 'bp','wl')
print m
print b
print r
We can also apply the barometric response function technique outlined by Rasmussen and Crawford (1997). This method uses convolution to least squares regression to identify and remove barometric effects.
In [9]:
def baro_eff(df,bp,wl,lag=100):
df.dropna(inplace=True)
dwl = df[wl].diff().values[1:-1]
dbp = df[bp].diff().values[1:-1]
#dwl = df[wl].values[1:-1]
#dbp = df[bp].values[1:-1]
df['j_dates'] = df.index.to_julian_date()
lag_time = df['j_dates'].diff().cumsum().values[1:-1]
df.drop('j_dates',axis=1,in_place=True)
# Calculate BP Response Function
## create lag matrix for regression
bpmat = tools.lagmat(dbp, lag, original='in')
## transpose matrix to determine required length
## run least squared regression
sqrd = np.linalg.lstsq(bpmat,dwl)
wlls = sqrd[0]
cumls = np.cumsum(wlls)
negcumls = [-1*cumls[i] for i in range(len(cumls))]
ymod = np.dot(bpmat,wlls)
## resid gives the residual of the bp
resid=[(dwl[i] - ymod[i]) for i in range(len(dwl))]
lag_trim = lag_time[0:len(cumls)]
return negcumls, cumls, ymod, resid, lag_time, dwl, dbp
In [10]:
negcumls, cumls, ymod, resid, lag_time, dwl, dbp = baro_eff(data, 'bp','wl',100)
plt.figure()
lag_trim = lag_time[0:len(negcumls)]
plt.scatter(lag_trim*24,negcumls, label='b.p. alone')
plt.xlabel('lag (hours)')
plt.ylabel('barometric response')
Out[10]:
In [11]:
plt.figure()
x = data.index.to_datetime()[1:]
y1 = resid
y2 = data['wl'][1:]
fig, ax1 = plt.subplots()
plt.title('Corrected Water Level')
ax2 = ax1.twinx()
ax1.plot(x,y1,c='r',label='Residual Water Level')
ax1.set_ylabel('Residual Water Level (ft)', color='r')
ax2.set_ylabel('Uncorrected Water Level (ft)', color='b')
ax2.plot(x,y2,c='b',label='Uncorrected Water Level')
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, loc=3)
Out[11]:
Clark, W.E., 1967, Computing the barometric efficiency of a well: Journal of the Hydraulics Division, American Society of Civil Engineers, v. 93, no. HY4, p. 93-98.
Merritt, M.L., 2004, Estimating hydraulic properties of the
floridan aquifer system by analysis of earth-tide, ocean-tide, and barometric effects, Collier and Hendry Counties, Florida: U.S. Geological Survey Water-Resources Investigations Report 03-4267, 70 p.
Rasmussen, T.C., and Crawford, L.A., 1997, Identifying and removing barometric pressure effects in confined and unconfined aquifers: Ground Water, v. 35, n. 3, pp. 502-511. doi: 10.1111/j.1745-6584.1997.tb00111.x